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Super-bubbles are shells in the interstellar medium produced by the simultaneous 
explosions of many supernova remnants. The solutions of the mathematical diffusion and 
of the Fourier expansion in ID, 2D and 3D were deduced in order to describe the diffusion 
of nucleons from such structures. The mean number of visits in the the case of the Levy 
flights in ID was computed with a Monte Carlo simulation. The diffusion of cosmic 
rays has its physical explanation in the relativistic Larmor gyro-radius which is energy 
dependent. The mathematical solution of the diffusion equation in ID with variable 
diffusion coefficient was computed. Variable diffusion coefficient means magnetic field 
variable with the altitude from the Galactic plane. The analytical solutions allow us to 
calibrate the code that describes the Monte Carlo diffusion. The maximum energy that 
can be extracted from the super-bubbles is deduced. The concentration of cosmic rays 
is a function of the distance from the nearest super-bubble and the selected energy. The 
interaction of the cosmic rays on the target material allows us to trace the theoretical 
map of the diffuse Galactic continuum gamma-rays. The streaming of the cosmic rays 
from the Gould Belts that contains the sun at it's internal was described by a Monte 
Carlo simulation. Ten new formulas are derived. 
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1. Introduction 

The spatial diffusion of cosmic rays (CR) plays a relevant role in astrophysics due 
to the fact that the binary collisions between CR and interstellar medium are neg- 
ligible, e.g. j n these last years a general consensus has been reached on the 
fact that particle acceleration in Supernova Remnants , in the following SNR, may 
provide CR at energies up to 10 15 eV^and probability density function in energy 
p(E)oc E~ 2 . Some of methods for the numerical computation of the propagation of 
primary and secondary nucleons have bee n imp lemented ^El The previous efforts 
have covered the diffusion of CR from 

snrEEI 

but the SNR are often concentrated 
in super-bubbles (SB) that may reach up to 800 pc in altitude from the Galactic 
plane. The complex 3D structure of the SB represents the spatial coordinates where 
the CR are injected. We remember that the super-bubbles are the likely source of at 
least a substantial fraction of GCRs; this can be shown from the analysis of the data 



f 
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from the Cosmic Ray Isotope Spectrometer (CRIS) aboard the Advanced Composi- 
tion Explorer (ACE) and in particular the 22Ne/20Ne ratio at the cosmic-ray source 
using the measured 21Ne, 19F, and 170 abundances as 'tracers" of secondary pro- 
duction^. Further on is known that w 75% of supernova occur in super-bubbles, 
and w 88% of the cosmic-ray heavy particles are acceler ated there because of the 
factor of 3 enhanced super-bubble core metallicity . 

In order to describe the propagation of CR from SB, in Section [2] we have 
developed the mathematical and Fourier solutions of the diffusion in presence of 
the stationary state and the Monte Carlo simulations which allow us to solve the 
equation of the diffusion from a numerical point of view . The case of random walk 
with variable length of the step was analysed in Section |2~B"1 The physical nature of 
the diffusion as well as evaluations on transit times and spectral index are discussed 
in Section [3] The 3D diffusion with a fixed length of a step and injection points 
randomly chosen on the surface of the expanding SB was treated in Section |4] 



2. Preliminaries 

In this section the basic equations of transport are summarised, the rules that govern 
the Monte Carlo diffusion are set up and the various solutions of the mathematical 
and Fourier diffusion are introduced. The statistics of the visits to the sites in the 
theory of Levy flights and it's connection with the regular random walk are analysed. 
It is important to point out that the mathematical and Monte Carlo diffusion here 
considered cover the stationary situation. 



2.1. Particle transport equations 

The diffusion-loss equation as deduced by^l ( equation (20.1) ) for light nuclei once 
the spallation phenomena are neglected has the form 

?J± = DV 2 N t + j^[b(E)N t ] + Q h (1) 

here Ni is the number density of nuclei of species i , V 2 is the Laplacian operator , 
D is the scalar diffusion coefficient, J^[b(E)Ni] takes account of the energy balance, 
and Qi represents the injection rate per unit volume. 

When only protons are considered and the energy dependence is neglected Equa- 
tion JTJ) becomes 

8P(x,y,z,t) 2 , . , , 
— = DV P(x,y,z,t) , (2) 

where P is the probability density , see^^ (equation 12.34b ). In this formulation 
the dependence for the mean square displacement R 2 (t) is, see^^ (equation. 8.38 ),: 



R 2 (t) = 2dDt (t -> oo) 



(3) 
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where d is the spatial dimension considered. From equation the diffusion coef- 
ficient is derived in the continuum : 

In the prese nce of discrete time steps on a 3D lattice the average square radius after 
N steps, see ^ (equation 12.5 ), is 

(R 2 (N))~2dDN , (5) 

from which the diffusion coefficient is derived 



If (R 2 (N)) ~ N , as in our lattice : 



when the physical units are 1 and 



D= 2^d XVtr ' (8) 

when the step length of the walker is A and the transport velocity v tr - 
Another useful law is Fick' s second equation in three dimensions ^El, 

f -DVO , (9) 

where the concentration C(x,y,z) is the number of particles per unit volume at the 
point (x,y,z). We now summarise how Fick' s second equation transforms itself when 
the steady state is considered in the light of the three fundamental symmetries. In 
the presence of spherical symmetry, equation ([9]) becomes 

1 d , o dC , „ / i 

r dr dr 

the circular symmetry (starting from the hollow cylinder) gives 

!fr§>-° • <"> 

and the point symmetry (diffusion from a plane) produces 



in the case of constant D and 



dr dr 



when D is a function of the distance or the concentration. 
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2.2. Monte Carlo diffusion 

The adopted rules of the d-dimensional (with d=l,2,3) random walk on a lattice 
characterised by NDIA4 d elements occupying a physical volume side d (in the hy- 
pothesis of a stationary state and fixed length step) are specified as follows: 

(1) The motion starts at the center of the lattice , 

(2) The particle reaches one of the 2 x d surroundings and the procedure repeats 
itself , 

(3) The motion terminates when one of the 2 x d boundaries is reached , 

(4) The number of visits is recorded on M. , a d-dimensional memory or concen- 
tration grid , 

(5) This procedure is repeated NTRIALS times starting from (1) with a different 
pattern , 

(6) For the sake of normalisation the d-dimensional memory or concentration grid 
M d is divided by NTRIALS . 

2.3. Mathematical diffusion 

The solutions of the mathematical diffusion are now derived in 3D, 2D, ID and ID 
with variable diffusion coefficient. 



2.3.1. 3D analytical solution 

Consider a spherical shell source of radius b between a spherical absorber of radius 
a and a spherical absorber of radius c, see Figure Q] that is adapted from Figure 3.1 

of El. 

The concentration rises from at r=a to a maximum value C m at r=b and then 
falls again to at r=c . The solution of equation (TIT))) is 

C(r)=A+- , (14) 
r 

where A and B are determined from the boundary conditions , 

C(r)=C m (l-^(l-^y 1 a<r<b , (15) 

and 

C(r) =C m (^-l) b<r<c . (16) 

These solutions can be found in ' ' or in 1- . The second solution (equation (fir?)) ) 
can be applied to explain the 3D diffusion from a central point once the following 
connection between continuum and discrete random walk is made 

(17) 



NDIM - 1 



s/,h NDIM + 1 



NDIM - 1 
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For example when NDIM ^> 1 or | ^> 1 

C(r)~C m - , (18) 
r 

in other words, the concentration scales as an hyperbola with unity represented 
by the mean free path b. The 3D theoretical solution as well as the Monte Carlo 
simulation are reported in Figure [21 

2.3.2. 2D analytical solution 

The 2D solution is the same as the hollow cylinder. The general solution to equa- 
tion (JTTJ) is 

C(r) = A + Bln{r) . (19) 




Fig. 1. The spherical source is represented by a dashed line, the two absorbing boundaries with 
a full line. 
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Fig. 2. Values of or concentration computed with equation l|16ll (full line) compared 

with the results of a Monte Carlo simulation (filled circles). The parameters are: NDIM=101, 
NTRIALS=100 and side=800pc. 



The boundary conditions give 

C{r) = cJ^\ a<r<b , (20) 
ln(b/a) 

and 

c ^= cJ ml "^^ ■ < 21 > 

The new 2D theoretical solution as well as the Monte Carlo simulation are reported 
in Figure [3l 



2.3.3. ID analytical solution 

The ID solution is the same as the diffusion through a plane sheet. The general 
solution to equation (fT2|) is 

C(r) = A + Br . (22) 

The boundary conditions give 

C(r) = C m 7 —^ a<r<b , (23) 
o — a 
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Fig. 3. Values of .M^ 2 ' or concentration computed with equation (12 Hi (full line) compared 
with the results of a Monte Carlo simulation (filled circles). The parameters are: NDIM=101, 
NTRIALS=4000 and side=800pc. 

and 

C{r)=C m r —^ b<r<c . (24) 
o — c 

The ID theoretical solution as well as the Monte Carlo simulation are reported in 
Figure 2] 



2.3.4. ID analytical solution with variable D 

The diffusion coefficient D is assumed to be a function ,f ,of the distance 

D = D {l + f(r)) , (25) 

where Dq is the value of D at the plane where the diffusion starts. On introducing 
the integral I defined as : 

(26) 



l + /(r) ' 

the solution to equation (fl~3|) is found to be , 

C-Ci h-I 



C1-C2 J1-/2 ' (2?) 
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Fig. 4. Values of .M^ 1 ' or concentration computed with equation i'2-li (full line) compared 
with the results of a Monte Carlo simulation (filled circles). The parameters are: NDIM=101, 
NTRIALS=100 and side=800pc. 



where C\ ,Ci , 1\ and I2 are the concentration and the integral respectively com- 



puted at r = r\ and r 
is assumed to be 



'"2 



For the sake of simplicity the diffusion coefficient 



D = D (l + a D (r- ri )) , (28) 

where ac is a coefficient that will be fixed in Section T3.2.2I . A practical way to 
determine arj is connected with the interval of existence of r , [r±, r 2 ], and with the 
value that D takes at the boundary 

D(r 2 )/D - 1 

a D = . (29) 

r 2 - ri 

Along a line perpendicular to the plane the concentration rises from at r=a to a 
maximum value C m at r=b and then falls again to at r=c . The variable diffusion 
coefficient is 



and 



D = D (l + a D {b-r)) 
when a < r < b , 

D = D (l + a D (r - b)) 
when b < r < c 



(30) 



(31) 
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Fig. 5. Values of jVf' 1 ' or concentration computed with equations l|32fl and l|33|l when a =- 
400 pc, b=0 pc and c=400 pc : D(c)/D(b) = 1.0001 ( full line ), D(c)/D(b) = 10 ( dashed ), 
D(c)/D(b) = 100 ( dot-dash-dot-dash ), D(c)/D(b) = 1000 ( dotted ). 



The solution is found to be 

_ (In (-a D a + 1 + a D b) - In (-a D r + 1 + a D bj) C m 

L/ y r > ~ 1 — 7 ; — ^ — ; r% \ 6Z ) 

m (— aD a + 1 + od b) 

when a < r < b , 

and 

C m (- In (cd c + 1 - a D b) + In (a D r + 1 - ao b)) . . 

C V ) = j— 7 r"i M ( 33 ) 

In (ao c + 1 — «d oj 

when b < r < c 

The new ID theoretical solution when the diffusion coefficient is variable is reported 
in Figure [5l 

Confront our Figure [5] with Figure 8 in^ , where the density of primary CR is a 
function of the distance from the Galactic plane and different physical parameters. 



2.4. Fourier Series 

In the case of restricted random walk the most straightforward way to compute 
the mean number of visits is via the solution of the Fokker-Planck equation, which 
in the c ase of simple, symmetric random walk equation reduces to the diffusion 
equation^^. Given the dimension d, the random walk takes place on a finite domain 
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D = [— L/2, L/2] d C R d , and the boundaries are supposed to be absorbing. In this 
case the solutions of the diffusion equation can be obtained with the usual methods 
of eigenfunction expansion an d they are cosine series whose coefficients decay 

exponentially, see equation (25) in^^l The solutions particularized as a cuts along 
one axis are 



with 



2 

M {1) (x) = - «m COS 

m=0 

M^{x,Q)=(\X £ k 



/COS 



(2m + 1)ttx 
L 

(2m + 1)ttx 



m,l=0 
3 oo 



^ K m lj COS 
m,l,j=0 



L 
2 



(2m + 1)ttx 
L 



(34) 
(35) 
(36) 



2L 2 



(2m + 1) 2 tt 2 



2L 2 



((m + l) 2 + (Z + l) 2 )7T 2 



(2L 2 ) 



((m + l) 2 + (/ + l) 2 + (j + l) 2 )7r 2 • 

These three new solutions are the Fourier counterpart of the solutions of the mathe- 
matical diffusion (|24|) , (|2"T|) and (fT6|) respectively. The comparison between solution 
in 3D with methods of eigenfunction expansion, mathematical diffusion and Monte 
Carlo simulation is carried out in Figure [S] 



2.5. Levy random walks in ID 

In the case of the Levy random walk in ID it should noted that 6x, i.e. the single 
step, can be written as 

Sx = £1, (37) 

where I > is the length of the flight and £ is a stochastic variable taking values 
1 and — 1 with equal probability. For the distribution of I chose 



p(l) = cr (Q+1) 



(38) 
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Fig. 6. Profile of M^ 3) versus the number of sites as given by Monte Carlo simulation with 
regular lengths step (filled circles ) + solution of the mathematical diffusion (full line ) + Fourier 
expansion (dot-dash-dot-dash line ) . The parameters are NDIM=101 , L=25, NTRIALS=10 000 
, m, I , j ranging from 1 to 21. 
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where < a < oo, and then 

p(6x) = l/2p(l) . (39) 

Let S a denote the average number of visits to a site obtained via a Monte Carlo 
simulation ; when Ai a denotes the corresponding values given by the equation 122, 

, AOll N 2 ^ L a . (rrm(x + a)\ . (mna\ 
M ^ = L^(^- a Sm {-^—) Sm {—) ' (40) 

n— 1 x ' N ' 

where 

D=-(-Yf Mr (41) 
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being 

L 

T = Y / S a (x) . (42) 



x=0 



Once we calculate the constant in pdf (I38|) (probability density function) the Pareto 

1991 

distribution, P, is obtained 1 ^ 

P(l;a,c) = -^ I , (43) 

where a < I < oo , a >0,c >0. The average value is 

I=^T , (44) 
c — 1 



^ 2 = „ — ^ . (45) 



which is defined for c > 1, and the variance is 

(-2 + c)(-l + c) 

which is defined for c > 2. In our case the steps are comprised between 1 and I 
and the following pdf , named Pareto truncated Pt , has been used 



fir(*;*mas,a)=( 1 )|^T ■ (46) 

The average value of Pt is 



- O; ( Ijnax H~ ^max) f AT\ 

= ~(l max a -!)(-! + a) ' 1 



and the variance of Pt is 

2 ™ 1 '"max \ ''max ^ ''max ™ ^ 



2a ±i a _ 9 / a + 2 ~, _ 9 / a + 1 rv 2 , 2/ a+2\ 

mas Lt i LX hnax I 



(-2 + a) (-1 + a) 2 (-/ ma:c 2a + 2 Z ma:c Q - l) 

rv f J a rv^ — 9 rv 1 Q _|_/l/ a ~^^rv / 2 _|_ 7 a + 2\ 

™ \ ''max ^ ^ ''max ~f ^ ''max ''max ~r hnax I , . _\ 

+ — 27 27i ~ s ■ ( 48 ) 

(-2 + a) (-1 + a) (-Uar " + 2/„ iQ2: Q -l) 

This variance is new formula and is always defined for every value of a > 0; con- 
versely the variance of the Pareto distribution can be defined only when a > 2. 
From a numerical point of view a discrete distribution of step's length is given by 
the following pdf 

p(l)dl oc r (Q+1) dl , (49) 

and the set of lengths could be obtained through a numerical computation of the 
inverse function <=Q. The interval of existence of the length-steps is comprised be- 
tween 1 (the minimum length) and Z maa; =NDIM that is number that characterises 
the lattice . We can now examine how the impact of the Levy random walk influ- 
ences the results of the simulations in ID; in the following NTRIALS will represent 
the number of different trajectories implemented in the Monte Carlo simulations. 
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pc 

Fig. 7. Profile of jVf' 1 ' versus the distance from the Galactic plane as given by Monte Carlo sim- 
ulation with regular lengths step + theoretical results given by equation i 1o4[ ) + Monte Carlo simu- 
lation with Levy random walk (filled circles , dotted line and dot-dash-dot-dash line respectively ). 
The parameters are NDIM=51 , L=25, m ranging from 1 to 20, side =800 pc , NTRIALS=10000 
and o=1.8 . 



The memory grid A^ 1 ) as a function of the distance from the center takes the cu- 
rious triangular shape visible in Figure [7] ; a theoretical explanation resides is the 
random walk in bounded domains EE3I. i n th a t. figure the theoretical solution ( for- 
mula (25) in when dimension d=l) is also reported and the fit is satisfactory. In 
Figure [7] we report the visitation grid as given by the presence of the Levy random 
walk once a certain value of a is chosen; is evident the effect of lowering the value of 
the visitation grid respect to the regular random walk. The effect of increasing the 
exponent a that regulates the lengths of the steps could be summarised in Figured] 
in which we report the value of the memory grid at the center of the grid. 

From the previous figure is also clear the transition from Brownian motion (a w 2 
) to regular motion (a s» 10 ). 



2.6. Levy diffusion coefficient in ID 

The ID diffusion coefficient in presence of Levy flights can be derived from analyt- 
ical or numerical arguments that originate four different definitions summarised in 
Figure [9] Let the single displacement, fe, be defined by 



5x = £Z, 



(50) 
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Fig. 8. Behaviour of Ai' 1 ' at the center of the lattice as a function of a (dashed— line) + theo- 
retical solution (dotted-line) + asymptotic behaviour at (NDIM+l)/2 (full line). Parameters as 
in Figure [7] . 



where I > is the length of the flight and £ is a discrete random variable charac- 
terised by the following probabilities p(l)=+l and p(2) =-1. The probability p{5x) 
of a single displacement Sx is 

p(8x) = ±p(l) . (51) 

This is the product of the probability = 1/2 to go at right or at left by the 
probability p(l) to make a step of length I. From (Sx) — and from the independence 
of the little steps follows 

(x(n)) = n(Sx) = , (52) 

where x(n) is the walker's position at the time n ( after n steps). From (Sx) = 
follows 

(S 2 x) =a 2 (Sx), (x 2 (n)) =a 2 (x(n)), (53) 
and due to the independence of the small steps , 

a 2 (x(n)) = na 2 (Sx) . (54) 
The connection with the diffusion coefficient is through the theoretical relationship 

D = i CT 2 (fa), (55) 
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Fig. 9. Behaviour of the diffusion coefficient as function a: DI. ( full line ), DU ( dashed ), 
D^ um ( dot-dash-dot-dash ), D^ l I urn ( dotted ). 



thus 



(56) 



This formula can be compared with the well know n fo rmula that can be found 
on the textbooks , see for example formula (12.5) in U^l We should now compute 
a 2 (5x). Since (Sx) = and £, I are random independent variables we have 

a*(Sx) = {6 2 x) = my) = {e)(l 2 ) ■ (57) 

Since (£ 2 ) = 1 

a 2 {8x) = (l 2 ) . (58) 

We can continue in two ways . A first one consists into evaluate I 2 according to the 
discrete distribution represented by equation (|49|) . This originates the first theoreti- 
cal diffusion coefficient , in the following D^ h . A second way consists into evaluate 
(I 2 ) in view of pdf (|46]) 



and consequently , 



(r 



DlL = 1/2 



(In 



l)(-2 + a) ' 



(/,, 



(U a - 1) (-2 + a) ' 



(59) 



(60) 
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where the new D*l eo is the second theoretical diffusion coefficient. 

The numerical techniques allow to implement two different definitions of the 
diffusion coefficient. The first one considers an infinite lattice and a variable number 
of Monte Carlo steps N (30, 31. ..35) each one connected with an index j. The diffusion 
coefficient in then computed as : 

v^6 R(i) 2 

j Li=l 2NW) (61) 

The second evaluation implements a typical Monte Carlo run on a lattice , see for 
example the data in Figure [7j Now R 2 is fixed and takes the value jv ~ p/ 9 M+1 in 
contrast with N , the number of iterations before to reach the two boundaries, that 
is different in each run. We therefore have 

DZ m = ^ , (62) 

where N is the average number of iterations necessary to reach the external bound- 
ary in each trial. Up to now in order to test the value of the diffusion coefficient 
we have fixed in one the minimum value that the step can take. We now explore 
the case in which the minimum length in pdf (|46p is i TO ,- n . On adopting the same 
technique that has lead to deduce Dlf leo we obtain a new more general expression , 
D th ' eo , for the ID diffusion coefficient 

n II,G _ i icy \ L min 'mai 'max 'min j 

U theo ~ ~ cW7 i a , ; ' \® 6 ) 



(o 2) ( /max ~l~ ^min ) 



or introducing the ratio r/ — 



2 (_„.2 



d u,g_ 1/2 °U (-n 2 + n") 
Uth ™ ( a -2)(-l+r,«) ■ 



3. Energy evaluations 

In order to introduce the physical mean free path we should remember that an 
important reference quantity is the relativistic ions gyro-radius, pz- Once the energy 
is expressed in 10 15 eV units ( E15) , and the magnetic field in 10 -6 Gauss ( H-q) 
we have: 

Pz = l.08-^-pc , (65) 

where Z is the atomic number. The solutions of the mathematical diffusion in pres- 
ence of the steady state do not require the concept of mean free path. The solutions 
of the diffusion equation trough a Fourier expansion, Levy flights and Monte Carlo 
methods conversely require that the mean free path should be specified like a frac- 
tion of the size of the box in which the phenomena is analysed . The presence of 
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irregularities allows the diffusion of the charged particles. The mean free path A sc , 
see for example 0^ , is 

A sc ~pz$- 2 , (66) 

with $ = -g^ where B\ are Bq are the strength of the random and mean magnetic 
field respectively. The type of adopted turbulence, weak or strong, the spectral index 
of the turbulence determine the exact correspondence between energy and mean free 
path. Some evaluation s on the relat ivistic mean free path of cosmic rays and Larmor 
radius can be found in I12I25I2£I2Z[ The transport of cosmic rays with the length of 
the step equal to the relativistic ion gyro-radius is called Bohm diffusion ^Hl and the 
diffusion coefficient will be energy dependent. The assumption of the Bohm diffusion 
allows to fix a one-to-one correspondence between energy and length of the step of 
the random walk. The presence of a short wavelength drift instability in the case 
of qua si-perpendicular shocks drives the diffusion coefficient to the "Bohm" limit 
29 30 rpjjg max i mum energies that can be extracted from SNR and Super-bubbles 
are now introduced. The transit times in the halo of our galaxy are computed in the 
case of a fixed and variable magnetic field. The modification of the canonical energy 
probability density function , p(E) cx E~ 2 , due to an energy dependent diffusion 
coefficient is outlined. 



3.1. Maximum available energies 

The maximum energy available in the accelerating regions is connected with their 
size, in the following A R and with the Bohm diffusion 

A R = Pz . (67) 

The most promising sites are the SNR as well as the super-shell. In the case of 
SNR, the analytical solution for the radius can be found in^J (equation 10.27) : 

R(i) =p^W) ,/5 , (68) 

where po is the density of the surrounding medium supposed to be constant , E expl 
is the energy of the explosion, and t$NR is the age of the SNR. Equation (|68|) can 
be expressed by adopting the astrophysical units 

pexpl ,2 

R = (^51_j4)l/5 12 . 47pc (6g) 

here £4 is £snr expressed in units of 10 4 years, Et^ pl is the energy expressed in u nits 
of 10 51 ergs and no is the density expressed in particles cm -3 . In agreement with^H 
Pa = no m was taken, where m = 1.4 win ■ this energy conserving phase end s at 
ti « 1.4. The thickness of the advancing layer is A R sa i?/12 according to^^ and 
this allows us to deduce the maximum energy ,£^ ax , that can be extracted from 
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SNR for a relativistic ion , see formula (f6"?| : 

E™ ax = 0.961 ( 51 4 ) /5 g_ 6 Z . (70) 
no 

This new formula is in agreement with more detailed evaluations E2l. 

The super-shells have been observed as expanding shells, or holes, in the HI- 

PT3l 

column density distribution of our galaxy - . The dimensions of these objects span 
from 100 pc to 1700 pc and often present elliptical shapes or elongated features, that 
can be explained as an expansion in a non-homogeneous medium . An explanation 
of the super-shell is obtained from the theory of the Superbubbles , in the following 
SB, whose radius is EI 



R 



25 



14tt 



1/5 /p D \ 1/5 
P 



t 3/5 , (71) 



where t is the considered time, i?sN is the rate of supernova explosions, and Eq 
the energy of each supernova; this equation should be considered valid only if the 
altitude of the OB associations from the Galactic ,zob=0, is zero, and only for 
propagation along the Galactic plane. 

When the astrophysical units are adopted the following is found 



R = 111.56 pc( 51 7 )*, (72) 



no 

where ty is the time expressed in 10 7 yr units, E^ xpl is the energy in 10 51 erg, n is 
the density expressed in particles cm -3 , and N* is the number of SN explosions in 
5.0 • 10 7 yr. A s an example, on inserting in equation (|72[) the typical parameters of 
GSH 238 E3 , EH pl = l , £ 7 =2.1 , n =l and N*=80.9 , the radius turns out to be 
419 pc in rough agreement with the average radius of GSH 238 , ~ 385 pc. 

The maximum energy ,E"^ ax ' sb , that can be extracted from a SB ( from for- 
mula l|67[)) in the Bohm diffusion is therefore: 

rpexpl at* 

£^«.rt = 8 . 6 l( fgi ) 1/5 t 3 7 /r °H_ 6 Z . (73) 
no 



Once the typical parameters of GSH 238 are inserted in equation ([73]) we obtain 
j^max.s _ 323.6 7J_g Z. The exact value of the magnetic field in the expanding 
layer of the SB is still subject of research , the interval i?_ 6 G {1, ...200} being 
generally accepted EH. The interested reader can compare our equation ([73]) with 
equation (15) in EH. 



3.2. Transit times 

The absence of the isotope 10 Be, which has a characteristic lifetime r r = 0.39 x 
10 7 yr , implies the following inequality on t, the mean cosmic ray residence time , 

T > T r . (74) 
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Given a box of side L with a proton inserted at the center, the classical time of 
crossing the box t c is : 

t c = ± . (75) 
Once the pc and yr units are chosen , c = 0.306— and 

t c = 1.63Zi yr , (76) 

where L\ is the side of the box expressed in lpc units. 

When the side of the box is Li=800, t c w 1304 yr, and t c <C r r , we should 
discard the diffusion of the cosmic rays through the ballistic model. 

We now analyse firstly the case in which the value of the diffusion coefficient is 
constant and secondly the case in which the diffusion coefficient is variable. 



3.2.1. Transit times when D is fixed 

From equations ([3|) and JHJ) we easily find the time , t, necessary to travel the 
distance R 

CfrCX 

where Vt r — <kr c . We continue inserting R = h , L = L\ -pc and t = t\ 10 7 yr , 

I? 

h = 0.815 x 10" 7 — i- . (78) 

CtrM 

Figure [10] reports a plot of the transit times versus the two main parameters Ai and 
Ctr together with the constrains represented by the lifetime of the 10 Be. In this 
case the lifetime of the phenomena ( cross-hatched region ) is tentatively fixed in 
h = 10. 

Conversely in presence of Levy flights (diffusion coefficient as given by for- 
mula [M]) the transit time t\ takes the form 



_ 7 L\ (a-2)(-l + n °) 
ctrXi a (ri 2 - ri a ) 



= -0.815 x HT 7 ^^ ^ -f^ ■ (79) 



The effect of the Levy flights is to increase the value of the transit times. Is inter- 
esting to realize that, in this astrophysical form of the Levy flights, the steps with 
transport velocity greater than the light velocity are avoided. 



3.2.2. Transit times when D is variable 

The value of the pressure in the Solar surroundings takes the value , 
Pi2l0~ 12 dyne cm" 2 I^SI w hen the pressure is expressed in 10" 12 dyne cm" 2 units 

Assuming equipartition between thermal pressure and magnetic pressure, the 
magnetic field Bq at the Galactic plane turns out to be 

B = V510" 6 pi2 Gauss . (80) 
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Fig. 10. Transit times expressed in 10 7 yr units as a function Ai when Li=800 : ctr = 0.2 ( 
full line), ctr = 0.4 ( dashed), ctr = 0.6 ( dot-dash-dot-dash), ctr = 0.8 ( dotted). The forbidden 
region is represented through cross-hatched lines. 



The behaviour of the magnetic field with the Galactic height z in pc is assumed to 
vary as 

B = J^- , (81) 
1 + a D z 

with Dq defined in Section f2. 3. 41 through the Bohm diffusion , 

In order to see how the presence of a variable diffusion coefficient modifies the transit 
times given by equation (|75|) we performed the following Monte Carlo simulation. 
The ID random walk has now the length of the step function of the distance from 
the origin. In order to mimic equation (|28p the step has length 

X = Xl{l + a D z) . (83) 

Typical behaviour of the influence of the variable D on transit times has been re- 
ported in Figure fTTI It is clear that an increasing D ( a decreasing B) lowers the tran- 
sit times; a rough evaluation predicts that transit times decrease cx D(A00pc) /D(0). 
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Fig. 11. Transit times expressed in 10 7 yr units as a function of D(400pc)/D(0) = 10 when 
Li=800, Ai=8 and ct r = 0.2. Assuming Z=l ,H-e=l we have ( Bohm diffusion) Ei5=7-38 



3.3. The spectral index 

At present the observed differential spectrum of cosmic rays scale as E~ 2J5 in the 
interval 10 10 eV - 5.0 10 15 eV and after 5 10 15 eV as £T 3 07 . 

Up to now the more interesting result on the prediction of the spec tral index 
resulting from particle acceleration in shocks , i.e. SNR or SB , is due tolSZEHHH 

The predicted differential energy spectrum of the high energy protons is 



N(E)dE oc E~ 2 dE 



(84) 



This is the spectral index where the protons are accelerated ; i.e. the region that 
has size ss j| , where R is the radius of the SNR or SB. Due to the spatial diffusion, 
the spectral index changes. The first assumption is that the transport occurs at 
lengths of the order of magnitude of the gyro-radius, the previously introduced 
Bohm diffusion. Two methods are now suggested for the modification of the original 
index as due to the diffusion : one deals with the concept of advancing segment/circle 
or sphere and the second one is connected with the asymptotic behaviour of the 3D 
concentration. 
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3.3.1. Index from the diffusion coefficient 

On considering a diffusion function of the dimensionality d (see formula (|SJl) , the 
square of the distances in presence of Bohm diffusion will be 

-2 0.3310 7 c tI .£i 5 ii 2 , OK . 
= -ff-^j pc ■ ( 85 ) 

Let us consider two energies -©15,1 and Ei$% with two corresponding numbers of 
protons (Z=l) , N± and N2. We have two expressions for the two advancing spheres 
(d=3) characterised by radius i?i(-Ei5,i) and ^(^15,2)- 

The ratio of the two densities per unit volume , A/2,1, gives the differential 
spectrum 

^,! = ^ = f (|^) 3/2 • (86) 

P2 iVl i/15,2 



But in the source 



and therefore 



when d=3 or 



f = (f^) 2 , (87) 



iV 2 .i oc (^i) 3 ' 5 oroci?- 3 - 5 , (88) 

-^45,2 



iv 2 ,i cx e-p+w , (89) 

when the concentration is function of d. This new demonstration simply deals with 
the concept of advancing segment/circle or sphere. 



3.3.2. Index from the 3D concentration 

Starting from the asymptotic behaviour of the 3D concentration , formula ()18j) the 
ratio of the two energy concentrations is possible to derive a new formula 

N 2 ,i oc E- {2+1) . (90) 



4. CR Diffusion in the SB environment 

Once equipped to calculate the spatial diffusion and temporal evolution of the CR 
attention now moves to the 3D evolution of CR once they have been accelerated in 
the advancing layer of super-shells. 



4.1. SB 

The super-shells have been observed as expanding shells, or holes, in the H I-column 
density distribution of our galaxy and the dimensions of these objects span from 
100 pc to 1700 pc 1 ^. The elongated shape of these structures is explained through 
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introducing theoretical objects n amed sup er-bubble (SB); these are created by me- 
chanical energy input from stars B0) 41 | 31 | 

The evolution of the SB's in the ISM can be followed once the following parame- 
ters are introduced t age , the age of the SB in yr, At, the time step to be inserted in 
the differential equations, £ burst j the time after which the bursting phenomena stops 
, N* , which is the number of SN explosions in 5.0 ■ 10 7 yr. An analytical solution 
of the advancing radius as a function of time can be found when the ISM has 
constant density. 

The basic equation that governs the evolution of the SB ^21 [ s momentum 
conservation applied to a pyramidal section, characterised by a solid angle, Afl,: 

±(AM j R j )=pI$ASl i , (91) 

where Af2j is the solid angle along a given direction and the mass is confined into 
a thin shell with mass AMj . The subscript j indicates that this is not a spherically 
symmetric system. 

In our case the density is given by the sum of three exponentials 

n(z) = ni e- z2 ' H ^ + n 2 e~ z2 / H * 2 + n 3 e^ Hs , (92) 

where z is the distance from the Galactic plane in pc, n\ — 0.395 particles cm" 3 , 
f/i=127 pc, 7i2=0.107 particles cm -3 , 772=318 pc, 77,3=0. 064 particles cm" 3 , and 
7^3=403 pc and this requires that the differential equations should be solved along 
a given number of directions each denoted by j . By varying the basic parameters, 
which are the bursting time and the lifetime of the source characteristic structures 
such as hourglass-shapes, vertical walls and V-shapes, can be obtained. The SB 
are then disposed on a spi ral structure as given by the percolation theory and the 
details can be found in 1 ^. 

The influence of the Galactic rotation on the results can be be obtained by 
introducing the law of the Galactic rotation as given by S3, 

y R (i? ) = 220(*l) - 382 kmscc- 1 , (93) 

here Rq is the radial distance from the center of the Galaxy expressed in pc. The 
original circular shape of the superbubble at z=0 transforms in an ellipse through 
the following transformation, T r , 

' xl = x + 0.264y t 
7 : { yt = y (94) 
zl = z, 

where y is expressed in pc and t in 10 7 yr units ^H. In the same way the effect of 
the shear velocity as function of the distance y from the center of the expansion, 
Vshiftiy) can be easily obtained on Taylor expanding equation (|93p 

V. Mf t(v)= 84.04-|- Km/s . (95) 
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This is the amount of the shear velocity function of y but directed toward x in a 
framework in which the velocity at the center of the expansion is zero and is new. 

4.2. Monte Carlo diffusion of CR from SB 

The diffusing algorithm here adopted is the 3D random walk from many injection 
points ( in the following IP), disposed on N$b super-bubbles. In this case a time 
dependent situation is explored ; this is because the diffusion from the sources 
cannot be greater than the lifetime of SB. The rules are : 

(1) The first of N SB SB is chosen. 

(2) The IP are generated on the surface of each SB. 

(3) The 3D random walk starts from each IP with a fixed length of a step. The 
number of visits is recorded on a 3D grid M. . 

(4) After a given number of iteration , ITMAX , the process restarts from (iii) 
taking into account another IP. 

(5) After the diffusion from the n SB , the process restarts from (i) and the (n+1) 
SB is considered. 

The spatial displacement of the 3D grid Ai(i,j,k) can be visualised through 
the ISO-density contours. In order to do so, the maximum value A4(i,j, k) max and 
the minimum value M.(i, j, k) min should be extracted from the three-dimensional 
grid. A value of the visitation or concentration grid can then be fixed by using the 
following equation: 

M(iJ,k) chosen = M(iJ,k) min + (M(i,j,k) max -M(i,j,k) min ) x coef , (96) 

where coef is a parameter comprised between and 1. This iso-surface rendering is 
reported in Fig[12] ; the Euler angles characterising the point of view of the observer 
are also reported . 

From the previous figure it is possible to see a progressive intersection of the 
contributions from different super-bubbles: once we decrease ITMAX , the spiral 
structure of spatial distribution of SB is more evident , see Figure [13] . 

The results can also be visualised through a 2D grid representing the concentra- 
tion of CR along the Galactic plane x-y , see Figure HU or ID cut , see Figure [T5l . 

From this figure the contribution from the SB belonging to the arms is clear 
; this profile should be a flat line when the mathematical diffusion is considered. 
Another interesting case is the cut along the Galactic height z , see Figure [16] from 
which the irregularities coming from the position of the SB on the arms are visible; 
this profile should have a triangular shape when the ID mathematical diffusion from 
the Galactic plane is considered. 
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Fig. 12. Three dimensional surface representing the surface brightness, front view. The ISO- 
density contours are fixed through coef=0.55. The number of super-bubbles Nsb is 260 and the 
time evolution is followed as in EH, Ai=172.94, IP=200 , ITMAX= 20000 , t D = 2.26 10 7 yr and 
ctr = 0.5. Assuming Z=l ,H—s=l we have ( Bohm diffusion) i?is=159.68 . 




4.3. The gamma emission from the Galactic plane 

The cosmic rays with an energy range 0.1 GeV < E < 410 5 GeV | 44 | 45 | can p roc j uce 
gamma ray emission ( 30AleV < E < 30GeT^) from the interaction with the target 
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The same as Figure [13] but ITMAX 



6.7 10 5 , and coef= 





material ' ' ' '. The gamma ray emissivity will therefore be proportional to the 
concentration of CR. 

A typical map of 7-emission from the Galactic plane is visualised in Figure [17] 
, where the additive property of the intensity of r adiation along the line of sight is 
applied. Using the algorithm of the nearest IP 43 it is also possible to visualise the 
7-emission in the Hammer- Aitof projection, see Figure [TBI 
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Fig. 14. Concentration of CR at the x-y (z=0) Galactic plane. Data as in Figure H"2l 

The observations with EGRET ^ provide cuts in the intensity of radiation that 
can be interpreted in the framework of the random walk. We therefore report in 
Figure [19] the behaviour of the intensity of the gamma ray emission proportional 
to concentration of CR when the distance from the Galactic plane is expressed in 
Galactic latitude; in the figure the stars represent the plot (30° > I > 10°) of SSI 



4.4. Monte Carlo diffusion of CR from the Gould Belt 

The physical parameters concerning the Gould Belt as deduced in SSI, are reported 
in Table [TJ The total energy is such as to produce results comparable with the 
observations and the kinematic age is the same as in SSI In order to obtain E to t = 
6. 10 51 erg with t hnl ' st = 0.015 10 7 yr, we have inserted N*= 2000. The time necessary 
to cross the Earth orbit , that lies 104 pc away from the Belt center, turns out to 
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Fig. 15. Cut along the Galactic plane , x-y, of the concentration of CR crossing the center. Data 
as in Figure [12] 



Size (pc 2 ) 466 • 746 at b=0 

Expansion velocity (km s ) 17 
Age(10 7 yr) 2.6 
Total energy (10 51 erg) 6 



be 0.078 10 7 yr that means 2.52 10 yr from now. The 2D cut at 2=0 of the 
superbubble can be visualised in Figure 1201 Our model gives a radial velocity at z 
=0 T4hoo=3.67km s _1 . The influence of the Galactic rotation on the direction and 
modulus of the field of originally radial velocity the transformation (|93|) is applied 
, see Figure 1211 A comparison should be done with Figure 5 and Figure 9 in ^1 

The results of the streaming of CR from the complex 3D advancing surface of 
the Gould Belt can be visualized through a 2D grid representing the concentration 
of CR along the Galactic plane x-y (z=0) , see Figurel22"l or through a perpendicular 
plane x-z (y=Q), see Figure l23l 

The variation of CR due to the presence of the Gould Belt are also analyzed in 
Section 3.2 in ^. 
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Fig. 16. Cut of the concentration of CR along a direction perpendicular to the Galactic plane 
and crossing the center. Data as in Figure [12] 

5. Conclusions 

We have studied the propagation of CR with the aid of mathematical diffusion in 
ID, 2D and 3D. Due to the nature of the problem , diffusion from the Galactic 
plane , the first results can be obtained in ID. The ID diffusion has also been 
analysed in the case of a variable diffusion coefficient: which in our astrophysical case 
consists in a magnetic field function of the altitude z from the Galactic plane. These 
analytical results allow to test the Monte Carlo diffusion. The physical explanation 
is the diffusion through the Larmor radius which means that the diffusion coefficient 
depends on energy. The calculation of the transit times of CR explains the absence 
of the isotope 10 Be. Further on the theoretical and Monte Carlo set up of the Levy 
flights allows to increase the applications of the random walk. For example from an 
observed profile of CR function of distance from the Galactic plane is possible to 
evaluate the concavity: the ID Levy flights predict profiles concave up and the ID 
regular random walk predicts no concavity. The observations with EGRET in the 
gamma range , see Figure 3 in ESI ; provide profiles that are concave up and therefore 
the Levy flights with the shape regulating parameter a can be an explanation of 
the observed profiles, see Figure [7l 

The CR diffusion in a super-bubble environment is then performed by evalu- 
ating the maximum energies that can be extracted (ss 10 11 eV) and performing a 
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Fig. 17. Intensity of the gamma-rays along the line of sight when the galaxy is face on. The 
parameters of the diffusion are Z=l , to = 0.6 10 7 yr , cj r = 0.5 , ii— 6=1 , Ai=9.22 

and ( Bohm diffusion) E = 8.5 10 6 GeV 




pc 



Monte Carlo simulation of the propagation of CR from SB disposed on a spiral-arm 
network. According to the numerical results, the concentration of CR depends on 
the distance from the nearest SB and from the chosen energy. 

This fact favours the local origin of CR SSI that comes out from an analysis of 
the results of KASCADE data^l From an astrophysical point of view, local origin 
means to model the diffusion from the nearest SNR or SB. Like a practical example 
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Fig. 18. Intensity of the gamma-rays along the line of sight in the Hammer-Aitof projection. 
The parameters are the same of Figure \T7\ 
Contouring on a galactic grid 




we now explore the case of the Gould Belt, the nearest SB , which has an average 
distance of 300 « pc from us. The concentration of CR at the source of this SB is 
, according to formula (fT8|) : 



C n 



C(300) 



300 pc 
Pz 



277 x C(300) 



E 



15 



(97) 



where C(300) is the measured concentration on Earth at the given energy. Is in- 
teresting to point out that the crossing time (25 My ago) of the expanding Gould 
Belt with the Earth with a consequent increases in the concentration of CR ca n 
explain one event of the fossil diversity cycle that has a periodicity of 62 My . 
The gamma-ray emissivity is different in the various regions of the galaxy "21 , and 
in order to explain this fact the theoretical map of gamma-ray from the Galac- 
tic plane has been produced. CR with energies greater tha t 10 18 eV^ , U ltra-High 
Energy (UHE) , are thought to be of extra-galactic origin p7 | 53 | 54 | 55 | 56 | 57 |) and 
therefore the diffusion from extra-galactic sources through the intergalactic medium 
is demanded to a future investigation. 
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Fig. 19. Theoretical gamma-ray emission as given by Monte Carlo simulation as function of the 
Galactic latitude when the cut is crossing the center and the observer is at 700 pc, data as in 
Figure [T71 The stars represent the observational cut. 
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Fig. 21. The stars represents the rotation-distorted section of the Gould Belt and the big star 
the Sun. The velocity field of the expansion modified by the shear velocity is mapped. The Galaxy 
direction of rotation is also shown. 
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Fig. 22. Concentration of CR diffusing away from the Gould Belt at the x-y (z=0) plane. Ai=5.51, 
IP=10000 , ITMAX= 200 , t D = 7204 yr and c tr = 0.5. Assuming Z=l ,H- 6 =1 we have ( Bohm 
diffusion) Bi 5 =5.08 . 
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Fig. 23. The same as Figure [23] but in the x-z (y=0) plane. 
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